#include <boost\geometry.hpp>
#include <boost\geometry\geometries\point.hpp>
#include <boost/geometry/geometries/box.hpp>
#include <boost/geometry/geometries/polygon.hpp>

#include <boost/geometry/index/rtree.hpp>


#include <cmath>
#include <vector>
#include <iostream>
#include <boost/foreach.hpp>
#include <boost/shared_ptr.hpp>

namespace bg = boost::geometry;
namespace bgi = boost::geometry::index;

void main2()
{
	typedef bg::model::point< double, 2, bg::cs::spherical_equatorial<boost::geometry::degree> > point;

	typedef bg::model::box<point> box;

	typedef bg::model::linestring<point> line;
	typedef boost::shared_ptr<line> shp;

	typedef std::pair<box, shp> value;

	// create the rtree using default constructor
	bgi::rtree< value, bgi::linear<16, 4> > rtree;
	//bgi::rtree< value, bgi::quadratic<16> > rtree;
	//bgi::rtree< value, bgi::rstar<16, 4> > rtree;

	std::cout << "filling index with polygons shared pointers:" << std::endl;

	// create some polygons and fill the spatial index
	for (int i = 0; i < 45; ++i)
	{
		for (int j = 0; j < 45; ++j)
		{
			shp p(new line());
			p->push_back(point((j + 0), (i + 1)));
			p->push_back(point((j + 1), (i + 0)));

			// calculate bounding box
			box b = bg::return_envelope<box>(*p);
			// insert new value
			rtree.insert(std::make_pair(b, p));
		}
	}

	// find values intersecting some area defined by a box
	box query_box(point(0.5, 0.5), point(1.5, 1.5));
	std::vector<value> result_s;

	for (int k = 0; k<100000; ++k) {
		result_s.clear();
		rtree.query(bgi::intersects(query_box), std::back_inserter(result_s));
	}

	// display results
	std::cout << "spatial query box:" << std::endl;
	std::cout << bg::wkt<box>(query_box) << std::endl;
	std::cout << "spatial query result:" << std::endl;
	BOOST_FOREACH(value const& v, result_s)
		std::cout << bg::wkt<line>(*v.second) << std::endl;

}